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Systems of particles interacting with “stealthy” pair potentials have been shown to possess in- 
hnitely degenerate disordered hyperuniform classical ground states with novel physical properties. 
Previous attempts to sample the infinitely degenerate ground states used energy minimization tech¬ 
niques, introducing algorithmic dependence that is artificial in nature. Recently, an ensemble theory 
of stealthy hyperuniform ground states was formulated to predict the structure and thermodynamics 
that was shown to be in excellent agreement with corresponding computer simulation results in the 
canonical ensemble (in the zero-temperature limit). In this paper, we provide details and justifi¬ 
cations of the simulation procedure, which involves performing molecular dynamics simulations at 
sufficiently low temperatures and minimizing the energy of the snapshots for both the high-density 
disordered regime, where the theory applies, as well as lower densities. We also use numerical sim¬ 
ulations to extend our study to the lower-density regime. We report results for the pair correlation 
functions, structure factors, and Voronoi cell statistics. In the high-density regime, we verify the the¬ 
oretical ansatz that stealthy disordered ground states behave like “pseudo” disordered equilibrium 
hard-sphere systems in Fourier space. The pair statistics obey certain exact integral conditions with 
very high accuracy. These results show that as the density decreases from the high-density limit, 
the disordered ground states in the canonical ensemble are characterized by an increasing degree 
of short-range order and eventually the system undergoes a phase transition to crystalline ground 
states. In the crystalline regime (low densities), there exist aperiodic structures that are part of the 
ground-state manifold, but yet are not entropically favored. We also provide numerical evidence 
suggesting that different forms of stealthy pair potentials produce the same ground-state ensemble 
in the zero-temperature limit. Our techniques may be applied to sample the zero-temperature limit 
of the canonical ensemble of other potentials with highly degenerate ground states. 


I. INTRODUCTION 

There has been long-standing interest in the phase 
behavior of many-particle systems in d-dimensional Eu¬ 
clidean spaces in which the particles interact with 
soft, bounded pair potentials [1-12]. Considerable at¬ 
tention has been devoted to the determination of the 
classical ground states (global energy minima) of such 
interactions [3, 6, 11, 12]. While typical interactions lead 
to unique classical ground states, certain special pair po¬ 
tentials are characterized by degenerate classical ground 
states—a phenomenon that has attracted recent atten¬ 
tion [12-22]. 

One family of such pair interactions are the “stealthy 
potentials” because their ground states correspond to 
configurations that completely suppress single scattering 
for a range of wave numbers. The Fourier transforms of 
these potentials are bounded and non-negative and have 
compact support [12], and hence they have correspond¬ 
ing direct-space potentials that are bounded and long 
ranged. Because of their special construction in Fourier 


* torquato@electron.princeton.edu 


space, hnding the ground states of stealthy potentials 
is equivalent to constraining the structure factor to be 
zero for wave vectors k contained within the support of 
the Fourier transformed potential [12], as will be sum¬ 
marized in Sec. II. In the case when the constrained 
wave vectors lie in the radial interval 0 < |k| < K, the 
stealthy ground states fall within the class of hyperuni¬ 
form states of matter [23] and can be tuned to have 
varying degrees of disorder. Disordered hyperuniform 
systems in general are of current interest because they 
are characterized by an anomalously large suppression 
of long-wavelength density fluctuations and can exist as 
equilibrium or nonequilibrium states, either classically or 
quantum mechanically [24-37]. Moreover, because disor¬ 
dered hyperuniform states of matter have characteristics 
that he between a crystal and a liquid [12], they are en¬ 
dowed with novel physical properties [18, 19, 38-46]. 

When a dimensionless parameter y, inversely propor¬ 
tional to the number density p and proportional to K'^ 
(size of the constrained region) is sufficiently small, the 
hyperuniform ground states are inhnitely degenerate and 
counterintuitively disordered (i.e., isotropic without any 
Bragg peaks) [12]. However, when y is large enough (p is 
sufficiently small), there is a phase transition to a regime 
in which the ground states are crystalline or highly or- 



2 


dered [13-15, 19]. For each spatial dimension d, there is 
a special value of X; Xmax) which the ground state is 
unique [47]. The unique ground state is the dual (recipro¬ 
cal lattice) of the densest Bravais lattice packing in each 
dimension [12]. In two and higher dimensions, as soon as 
X drops below Xmaxi the set of the ground states become 
uncountably inhnite and gradually includes progressively 
less ordered structures [12]. Similarly to stealthy poten¬ 
tials, a family of two-, three-, and four-body potentials 
that lead to disordered ground states has also been de¬ 
fined in Fourier space and studied [17, 18, 21]. 

Due to the complexity of the problem, almost all pre¬ 
vious investigations of the ground states employed com¬ 
puter simulations. Such numerical studies were carried 
out in one, two and three dimensions [13, 14, 17-19]. The 
ground states were sampled by minimization of poten¬ 
tial energy at fixed densities starting from random initial 
conditions in a d-dimensional cubic simulation box under 
periodic boundary conditions. A few optimization tech¬ 
niques were employed to find the global energy minima 
with very high precision [14, 17]. 

Generally, a numerically obtained ground-state conhg- 
uration depends on the number of particles N within the 
fundamental cell, initial particle configuration, shape of 
the fundamental cell, and particular optimization tech¬ 
nique used [12]. Adding to the complexity of the problem 
is that the disordered ground states are highly degener¬ 
ate with a configurational dimensionality that depends 
on the density, and there are an infinite number of dis¬ 
tinct ways to sample this complex ground-state manifold, 
each with its own probability measure. These nontrivial 
aspects had made the task of formulating a statistical- 
mechanical theory of stealthy degenerate ground states 
a daunting one. Recently, we have formulated such an 
ensemble theory that yields analytical predictions of the 
structural characteristics and other properties of stealthy 
degenerate ground states [12]. A number of exact results 
for the thermodynamic and structural properties of these 
ground states were derived that applied to general en¬ 
sembles. We then specialized our results to the canonical 
ensemble (in the zero-temperature limit) by exploiting an 
ansatz that stealthy disordered ground states (for suffi¬ 
ciently small x) behave remarkably like “pseudo” disor¬ 
dered equilibrium hard-sphere systems in Fourier space. 
Our theoretical predictions for the pair correlation func¬ 
tion g 2 (r) and structure factor S{k) of these entropi- 
cally favored disordered ground states were shown to be 
in agreement with corresponding computer simulations 
across the first three space dimensions. We also made 
predictions for the corresponding excited states for suffi¬ 
ciently small temperatures that were in agreement with 
simulations. 

Because the focus of that previous investigation was 
the development of ensemble theories, few simulation 
details were presented about how the canonical ensem¬ 
ble was sampled to produce stealthy disordered ground 
states. One aim of the present paper is to provide a com¬ 
prehensive description of the numerical procedure that 


we used to produce the simulation results in Ref. 12. 
Moreover, here we also extend those results by apply¬ 
ing the simulation procedure to study numerically the 
ground states in the canonical ensemble for all allowable 
values of x and thus investigate the entire phase dia¬ 
gram for the entropically favored states across the first 
three space dimensions. In the second paper of this se¬ 
ries, we will study the exotic aperiodic “wavy phases” 
identihed in previous numerical work [14] (or “stacked- 
slider phases,” as called in the sequel to this paper [48]), 
a special part of the ground-state manifold. An analyti¬ 
cal model will enable an even more detailed study of this 
phase. 

As a justification of sampling the canonical ensemble 
instead of minimizing energy, we also demonstrate here 
how a variety of different optimization techniques affect 
the ground states that are sampled, which was not previ¬ 
ously investigated [14, 17, 18]. This investigation reveals 
that the pair statistics of the ground-state configurations 
indeed generally depend on the algorithm. Moreover, we 
show here that the energy minimization results depend 
on the initial conditions as well. We also provide the rea¬ 
son why the simulations in Ref. 12 and this paper employ 
noncubic, possibly deforming, simulation boxes for d > 2. 
Because almost all previous numerical simulations were 
performed using some specific form of stealthy potentials, 
we show here that different forms of stealthy potentials 
produce identical pair correlation functions, suggesting 
that the specihc choice of the potential form does not 
affect the ensemble being sampled. 

Among our major findings, we show that energy min¬ 
imizations starting from random initial conditions may 
lead to clustering of particles, the degree of which de¬ 
pends on the algorithm for a finite range of x below 1 /2 
across the hrst three space dimensions. When minimiz¬ 
ing the energy starting from configurations equilibrated 
at some temperature Te, the ground-state configurations 
discovered depend on Te- However, the algorithm depen¬ 
dence diminishes in the Te ^ 0 limit. We also demon¬ 
strate that the pair statistics [g 2 ir) and S{k)] in this limit 
do not depend on the particular form of the stealthy po¬ 
tential. The similarity between the structure factor in 
this limit and the pair correlation function of an equilib¬ 
rium hard-sphere system in direct space [12] is valid for 
X up to some dimension-dependent values between 0.25 
and 0.33 in the hrst three space dimensions. Beyond 
this range of Xj the hard-sphere analogy in Fourier space 
undergoes modihcation. As x increases further (to the 
value of about 0.4 in two dimensions, for example), the 
hrst peak in the structure factor diminishes while second 
peak in the structure factor grows and engulfs the hrst 
peak. Our simulated pair statistics obey certain exact 
integral conditions in Ref. 12 with very high accuracy, in¬ 
dicating the high hdelity of the numerical results. In the 
inhnite-system-size limit, at x = 0.5, the entropically fa¬ 
vored ground states undergo a transition from disordered 
states to crystalline states. Depending on the dimension, 
this phase transition can occur when aperiodic structures 
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still are part of the ground state manifold, demonstrat¬ 
ing that crystalline (ordered) structures can have a higher 
entropy than disordered structures. 

The rest of the paper is organized as follows: In Sec. II, 
we briefly summarize the numerical collective-coordinate 
procedure and other details of the simulation that we em¬ 
ploy in the present paper with justifications. In Sec. Ill, 
we study the dependence of the results on a variety of en¬ 
ergy minimization algorithms, initial conditions, and the 
forms of the stealthy potentials. In Sec. IV, we provide 
pair correlation function, structure factor, Voronoi cell- 
volume distribution, and configuration snapshots of the 
stealthy hyperuniform ground states obtained from the 
canonical ensemble in the zero-temperature limit. We 
provide concluding remarks and discussion in Sec. V, in¬ 
cluding suggestions for sampling the canonical ensemble 
in the zero-temperature limit of other potentials with de¬ 
generate disordered ground states. 

II. MATHEMATICAL RELATIONS AND 
SIMULATION PROCEDURE 

As detailed in Sec. II of Ref. 12, we simulate point 
processes in periodic fundamental cells (i.e. simulation 
boxes) with a pairwise additive potential v(r) such that 
its Fourier transform exists. Under nearest image con¬ 
vention, the total potential energy can be calculated by 
summing over all pairs of particles: 

( 1 ) 

i<3 

where N is the number of particles, = ri, r 2 ,..., r^v is 
the locations of the particles in d-dimensional Euclidean 
space, and = r—r^ . Instead of summing over all pair¬ 
wise contributions in the real space, the potential energy 
can also be represented in Fourier space: 


In a finite-sized system, the real-space pair potential has 
the same periodicity as the fundamental cell. Therefore, 
in the infinite-volume limit, the cell periodicity disap¬ 
pears. 

A family of “stealthy” potentials, which completely 
suppress single scattering for all wave vectors within 
a specihc cutoff in their ground states, are defined as 
[13, 14, 17-20]: 

U(fc), if|k|<A:, 

0 , otherwise, 

where V(k) is a positive isotropic function and AT is a 
constant. In this paper we always take K = 1, which 
sets the length scale. We will also use V{k) = 1 unless 
otherwise specihed. In the infinite-system-size limit, the 
isotropic f;(k) correspond to an isotropic real-space pair 
potential u(r) [12]. However, for finite systems, the cor¬ 
responding v(r) is anisotropic. In Appendix A, we com¬ 
pare the infinite-system-size limit u(r) with the hnite-size 
u(r)’s in different-shaped simulation boxes and select the 
simulation box shape to be used in this paper based on 
which u(r) is closest to the infinite-size-limit v(r). 

From Eqs. (2) and (5), one can see that a configuration 
is a stealthy ground state if h(k) = 0 for all k points such 
that 0 < |k| < A. Therefore, finding a ground state of a 
stealthy potential is equivalent to constraining h(k) = 0 
for all of those k points. However, in a simulation, one 
does not need to check all of the constraints. As detailed 
in Ref. 12, if there are (2M -|- 1) k points within the 
constrained radius, only M of them are independent and 
needed to be constrained to zero. Equation (2) can be 
simplified as [49]: 

<i>(r^) = —^u(k)Kk)|2 + $o, (6) 

where the sum is over all independent constraints, and 


$(r^) 


1 

2vf 


^{i(k)|h(k)|2 


iV^{i(k) , 

k 


( 2 ) 


where vp is the volume of the fundamental cell, u(k) = 
v(r) exp(—fk • r)dr is the Fourier transform of the 

pair potential, h(k) = exp(—zk • r^) is the complex 

collective density variable [with h(k = 0) = N], and 
both summations are over all reciprocal lattice vector k’s 
appropriate to the fundamental cell. For every k 7 ^ 0 , 
h(k) is related to the structure factor, S'(k), via 


S{k) 


|n(k)p 

N 


( 3 ) 


Given a v(k), the corresponding real-space pair potential 
is 


v{r) 


— {;(k) exp(zk 


r). 


( 4 ) 


$0 = [Af(iV-l)-2A^7(k)]/(2uF) (7) 

k 


is a constant independent of the particle positions . 
We now introduce a parameter 


M 

d{N-l)' 


( 8 ) 


which determines the degree to which the ground states 
are constrained, and therefore the degeneracy and disor¬ 
der of the ground states [14]. Note that the constraints 
depend on K and the fundamental cell but are indepen¬ 
dent of the specific shape of 7(k) as long as 7(k) > 0 
for all 0 < |k| < K. Therefore, changing {;(k) does not 
change the set of the ground states. However, there is 
no proof that changing ^(k) does not change the relative 
sampling weights of the ground states. 

In this paper we study various systems with different 
x’s and N's. One numerical complication is that these 
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numbers cannot be chosen arbitrarily, since M = xd{N — 
1 ) must be an integer consistent with the specific shape 
of the simulation box. (For example, a list of the allowed 
M values for a two-dimensional square box is given in 
Table II of Ref. 14.) This constraint is especially hard 
to meet when simulating multiple systems at the same x 
value across dimensions. In fact, both x and N in Table I 
(see Appendix C) had to be chosen carefully to meet this 
constraint. 

Taking the gradient of Eq. (6) yields the forces on par¬ 
ticles: 

F, = kf;(k) Im[n(k) exp(zk • r,)], 

k 

. 

where the sum is also over all independent constraints. 
This equation enables us to perform both energy mini¬ 
mizations and molecular dynamics (MD) simulations. In 
an energy minimization, a derivative-based algorithm is 
used. The first term on the right side of Eq. (6) is pro¬ 
vided to the algorithm as the objective function and the 
negative of the force in Eq. (9) is provided as the deriva¬ 
tive. In order to minimize energy, we have tried dif¬ 
ferent algorithms including the MINOP algorithm [50], 
the steepest descent algorithm allowing large steps [51], 
the low-storage BEGS (L-BFGS) algorithm [52-54], the 
Polak-Ribiere conjugate gradient algorithm [51, 55], and 
our “local gradient descent” algorithm described in Ap¬ 
pendix B. When x < 0-5, the objective function always 
ends up being very close to zero (the minimum). The 
maximum ending objective function for different algo¬ 
rithms varies from as high as 10“^ for a conjugate gra¬ 
dient algorithm to 10“^^ for the local gradient descent 
and steepest descent algorithms to 10“^° for the L-BEGS 
algorithm, and to as low as 10“^^ for the MINOP algo¬ 
rithm. Prom our practical point of view, all of these 
algorithms are precise enough, since an error of 10“^ or 
lower is indiscernible from any results presented below. 
Because the L-BEGS algorithm is the fastest, we will use 
it unless otherwise specified. 

The energy minimizations, if started from random ini¬ 
tial configurations, will sample an algorithm-dependent, 
nonequilibrium ensemble. To sample the canonical en¬ 
semble at a given equilibrium temperature Te we use 
MD simulations. One important parameter in MD sim¬ 
ulations is the integration time step. Since the optimal 
choice of the time step depends on the temperature, and 
the latter varies across several orders of magnitude in 
this paper, we desire a systematic way to determine the 
optimal time step. Starting from an energy minimized 
configuration and a very small time step (0.01 in dimen¬ 
sionless units), we repeat the following steps 10^ times to 
equilibrate the system and find a suitable time step: 

• Assign a random velocity from Boltzmann distri¬ 
bution at Te to each particle. 

• Calculate the total (kinetic and potential) energy 
of the system Ei. 


• Evolve the system 1500 time steps using the veloc¬ 
ity Verlet algorithm [56]. 

• Calculate the total energy of the system E^. 

• If I In ^ I > lx 10“^, then the time step is too 
large and errors will build up quickly. Therefore, 
we decrease the time step by 5%. On the other 
hand, if | In f I <4x 10 there is still some room 
to increase the time step. Since increasing the time 
step increases the efficiency of MD simulations, we 
increase the time step by 5%. 

After the system is equilibrated and the time step is cho¬ 
sen, we perform constant temperature MD simulations 
with particle velocity resetting [57]. A randomly cho¬ 
sen particle is assigned a random velocity, drawn from 
Maxwell-Boltzmann distribution, every 100 steps. We 
take a sample configuration every 3000 time steps until 
we have sampled 20 000 configurations unless otherwise 
specified. This amounts to an implementation of the gen¬ 
eration of configurations in the canonical ensemble. 

The above MD procedure works well for y < 0.5. How¬ 
ever, two new features arise when it is applied to x ^ 0-5 
in all dimensions. First, the potential energy surface de¬ 
velops local minima and energy barriers that can trap the 
system if Te is too small. We address this problem by 
using simulated annealing, employing a thermodynamic 
cooling schedule [58] which starts at T = 2 x 10“^ and 
ends at 10“®. Note that, by adopting a cooling schedule, 
we concede that we may only take one sample at the end 
of each MD trajectory, whereas a fixed-temperature MD 
trajectory produces multiple samples. 

The second new feature is that the entropically favored 
ground states are crystalline for x > 0-5. Unlike disor¬ 
dered structures, a crystalline structure has long-range 
order and may not “fit” in simulation boxes with certain 
shapes. To overcome the second problem, we simulate an 
isothermal-isobaric ensemble with a deformable simula¬ 
tion box. Every 20 MD time steps, 10 Monte Carlo trial 
moves to deform the simulation box are attempted. The 
pressure is calculated from Eq. (41) of Ref. 12. 

We employed the Wang-Landau Monte Carlo [59] to 
attempt to determine the entropically favored ground 
states for x > 0.5 in two and three dimensions. The 
Wang-Landau Monte Carlo is used to calculate the mi- 
crocanonical entropy §(<!>) as a function of the potential 
energy <I>. We limit our simulations to the energy range 
3 X 10“^° < 4) — <I>o < 10“® (in dimensionless units), 
where $0 is the ground state energy, by rejecting any 
trial move that violates this energy tolerance. This en¬ 
ergy range is evenly divided into 1000 bins. Starting from 
a perfect crystal structure in a simulation box shaped like 
a fundamental cell, small perturbations are introduced so 
the energy is within the range. After that, 60 stages of 
Monte Carlo simulations are performed, each stage con¬ 
taining 3 X 10^ trial moves. The “modification factor” in 
Ref. [59] is / = exp[5/(n -|- 10)], where n is the number 
of stages. 
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III. DEPENDENCE ON ENERGY 
MINIMIZATION ALGORITHM, MD 
TEMPERATURE, AND v{k) 

In this section, we present numerical simulation results 
demonstrating that: 


• Energy minimizations starting from Poisson initial 
configurations using different algorithms can yield 
ground states with different pair correlation func¬ 
tions. 


• Energy minimizations starting from MD snapshots 
at different temperatures can yield ground states 
with different pair correlation functions. 


• For configurations obtained by minimizing energy 
starting from MD snapshots at sufficiently small 
temperature, pair correlation functions do not de¬ 
pend on the minimization algorithm and the form 
of the stealthy potential. 


These results motivate the reason why we ultimately 
study and report results in Sec. IV in the canonical en¬ 
semble in the zero-temperature limit. For concreteness 
and visual clarity, we present results here in two dimen¬ 
sions. However, we have verified that all of the conclu¬ 
sions here also apply to one and three dimensions. 

We performed energy minimizations starting from 
Poisson initial configurations (i.e., Tg —)■ oo state at fixed 
density) using each of the five numerical algorithms men¬ 
tioned in Sec. II at y = 0.2 and y = 0.4. The results are 
shown in Figs. 1 and 2. At % = 0.2, the pair correla¬ 
tion functions produced by the MINOP algorithm and 
the L-BFGS algorithm are almost identical. However, 
the pair correlation function produced by the conjugate 
gradient algorithm noticeably differs. The steepest de¬ 
scent algorithm and our local gradient descent algorithm 
produce a significantly different pair correlation function 
with a much weaker peak at r = 0. The pair correla¬ 
tion functions produced by some algorithms appear to 
have g 2 {r) oc log(r) divergence near the origin. Since 
this divergence means particles have a tendency to form 
clusters, we call it a “clustering effect.” At y = 0.4, the 
clustering effect disappears, but the pair statistics pro¬ 
duced by different algorithms still differs. The fact that 
different optimization algorithms produce different pair 
statistics means that they sample the ground-state mani¬ 
fold with different weights. In other words, different opti¬ 
mization algorithms are sampling different ground-state 
ensembles. 



FIG. 1. (Color online) Pair correlation function as obtained 
from different optimization algorithms (as described in the 
legend) starting from Poisson initial configurations in two di¬ 
mensions at y = 0.2. Each curve is averaged over 20 000 
configurations of 136 particles each. The left inset zooms in 
near the origin, showing the differences between the five al¬ 
gorithms more clearly. The right inset uses a semilogarithmic 
scale to show 52 (r) oc log(r) near the origin. 



FIG. 2. (Color online) As in Fig. 1, except that y = 0.4 
and each curve is averaged over 20 000 configurations of 151 
particles each. The inset zooms in near the first well, showing 
the differences between the five algorithms more clearly. 

In order to avoid the complexity caused by the details 
of various optimization algorithms, we turn our interest 
to the canonical ensemble in the T —>■ 0 limit. To sam¬ 
ple this ensemble, we perform MD simulations at suffi¬ 
ciently small temperature Tgj, periodically take “snap¬ 
shots,” and then use a minimization algorithm to bring 
each snapshot to a ground state. To determine a “suf¬ 
ficiently small” Te, we calculated the pair correlation 
functions at various Te’s and present them in Fig. 3. The 
energy minimization result starting from Te —t 00 initial 
configurations clearly display the “clustering effect” at 
y = 0.2. When Te goes to zero, the “clustering effect” 
also diminishes. At y = 0.4, particles develop hard cores 
































6 


[52(0) = 0], therefore there is no clustering even if Te is 
large or infinite. However, the peak height of 52 (^) be¬ 
comes dependent on Te at this x value. For both x val¬ 
ues, the pair correlation functions of the two lowest Tj;’s 
are almost identical, verifying that the Te —>■ 0 limit ex¬ 
ists. These results show that Te = 2 x 10“® is sufficiently 
small in two dimensions. Similarly, we have found that 
Te = 2 X 10“^ and Te = 1 x 10“® are sufficiently small 
in one and three dimensions, respectively. These temper¬ 
atures are used in generating all of the results presented 
in Sec. IV A. 



(a) 



(b) 


FIG. 3. (Color online) Pair correlation function produced by 
L-BFGS algorithm starting from snapshots of MD at different 
equilibration temperatures Te, (a) x = 0.2 and (b) x = 0.4. 
Each curve is averaged over 20 000 configurations of 136 par¬ 
ticles each or 151 particles each. 

The energy minimization result starting from Poisson 
initial conhgurations differs for different algorithms, but 
the canonical ensemble in the T —>■ 0 limit should not 
depend on any particular algorithm. After finding that 
Te = 2 X 10“® is sufficiently small, we confirm the dis¬ 
appearing of algorithmic dependence by calculating the 
pair correlation function produced by different energy 
minimization algorithms starting from MD snapshots at 


Te = 2 X 10 Figure 4 shows the results. The curves 
for all algorithms almost coincide. 



FIG. 4. (Color online) Pair correlation function produced by 
the five different algorithms starting from snapshots of MD 
at equilibration temperature Te = 2 x 10“® at x = 0.2. Each 
curve is averaged over 20 000 configurations of 136 particles 
each. 


Last, the function V{k) in Eq. (5) can have different 
forms. This paper mainly use V(k) = 1 but we also 
want to know if the results obtained using this form 
are equivalent to those generated using other positive 
isotropic forms of V(k) as well. In principle, stealthy po¬ 
tentials of any form should have the same set of ground- 
state configurations, but the form of the stealthy poten¬ 
tial could theoretically affect the curvature of the poten¬ 
tial energy surface near each ground-state configurations 
and thus also affect their relative weights. Figure (5) 
shows the pair correlation function produced by different 
V{k)^s. The pair correlation functions for V{k) = 1 and 
V{k) = (I — k)^ at Te = 2 x 10“® are almost identical. 
For V{k) = (I — fc)®, we initially tried Te = 2 x 10“® 
but found that the “clustering effect” is still noticeable. 
We further lowered the temperature to Te = 2 x I0“^° 
to completely suppress the “clustering effect” to produce 
a pair correlation function identical to that of V(k) — I 
and V{k) = (I —fc)^ potentials. This result suggests that 
the functional form of V(k) does not produce noticeable 
differences in the ground-state ensembles in the T —)■ 0 
limit of the canonical ensemble. 
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FIG. 5. (Color online) Pair correlation function produced by 
different potentials starting from snapshots of MD at suffi¬ 
ciently low temperature at x = 0.2. Each curve is averaged 
over 20 000 configurations of 136 particles each. 


IV. CANONICAL ENSEMBLE IN THE T ^ 0 
LIMIT 

We will show here that the entropically favored ground 
states in the canonical ensemble in the T —>■ 0 limit for 
the first three space dimensions differ markedly below 
and above x = 0.5. For y < 0.5, the entropically favored 
ground states are disordered while for y > 0.5 the entrop¬ 
ically favored ground states are crystalline. Therefore, we 
will characterize them differently. For x < 0.5, we will 
report the pair correlation function, structure factor, and 
Voronoi cell statistics. For sufficiently small x, we will 
show that the simulation results agree well with theory 
[12]. For X > 0.5, we will report the crystal structures. 
The numbers of particles in all of the systems reported 
in this section are collected in Appendix C. 

A. X < 0-5 region 

Representative entropically favored stealthy ground 
states in the first three space dimensions at x = 0.1 and 
X = 0.4 are shown in Figs. 6-8. As x increases from 0.1 to 
0.4, the stealthiness increases, accompanied with a visu¬ 
ally perceptible increase in short-range order. This trend 
in short-range order is consistent with previous studies 
[14, 17, 18]. 



FIG. 6. (Golor online) Representative one-dimensional entropically favored stealthy ground states at (a) x = 0-1 S'^d (b) 
X = 0.4. 



FIG. 7. (Color online) Representative two-dimensional entropically favored stealthy ground states at (a) x = 0.1 and (b) 
X = 0.4. 



















FIG. 8. (Color online) Representative three-dimensional entropically favored stealthy ground states at (a) x = 0.1 and (b) 
X = 0.4. 



FIG. 9. (Color online) Structure factors for 1 < d < 3 for 0.05 < X ^ 0.33 from simulations and theory [12]. The smaller x 
simulation results are also compared with the theoretical results in the infinite-volume limit [12]. For x < 0.1, the theoretical 
and simulation curves are almost indistinguishable, and the structure factor is almost independent of the space dimension. 
However, simulated S{k) in different dimensions become very different at larger x- Theoretical results for x > 0.25 are not 
presented because they are not valid in this regime. 
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FIG. 10. (Color online) Pair correlation functions for 1 < d < 3 for 0.05 < X ^ 0-33 from simulations and theory [12]. The 
smaller x simulation results are also compared with the theoretical results in the infinite-volume limit [12]. For x < 0.1, the 
theoretical and simulation curves are almost indistinguishable. Theoretical results for x > 0.25 are not presented because they 
are not valid in this regime. 


We have calculated the pair correlation functions and 
the structure factors for various x values. Results for 
0.05 < X ^ 0-33 are shown in Figs. 9 and 10. The x < 0-2 
results are in excellent agreement with the “pseudo-hard- 
sphere ansatz,” which states that the structure factor 
behaves like pseudo equilibrium hard-sphere systems in 
Fourier space [12]. However, the theory gradually be¬ 
comes invalid as x increases. 

The pair correlation functions of the entropically fa¬ 
vored stealthy ground states are shown in Fig. 10. When 
X < 0.2, since the structure factor is similar to the pair 
correlation function of the hard-sphere system, inversely 
the pair correlation function is also similar to the struc¬ 
ture factor of the hard-sphere system. As x grows larger, 
the pseudo hard-sphere ansatz gradually deviates from 
the simulation result. 

We have checked that these pair statistics are consis¬ 
tent with four theoretical integral conditions of the pair 
statistics in the infinite-volume limit [12]. The first three 
conditions are Eqs. (58), (59), and (63) of Ref. 12, which 
are 

[ P{r)dr = 0, (10) 

f P{r)v{r)dr = 0, (11) 

./R'i 

and 

52 ( 0 ) = 1 - 2(ix + 2d^X / k'^~^Q{k)dk, (12) 

Jk 

where P{r) is the inverse Fourier transform of 0(fc — 
l)Q(fc), Q{x) is the Heaviside step function, and Q{k) = 
S{k) - 1. 


The fourth condition is that the pressure calculated 
from the “virial equation” [12] has to be either noncon- 
vergent or convergent to the pressure calculated from the 
energy route [12]. All pair statistics in Figs. 9 and 10 were 
generated using the step-function potential [the V{k) = 1 
case of Eq. (5)], but this potential does not lead to a con¬ 
vergent virial pressure. However, as we have shown ear¬ 
lier, the stealthy ground states that we generated here 
are also the ground states of other stealthy functional 
forms t'(k). In one dimension, to test our simulation 
procedure, we used the potential form V{k) = (1 — fc) 
to calculate the pressure from both the virial equation 
(Eq. (43) of Ref. 12) and the energy equation (Eq. (41) 
of Ref. 12). The pressure from the virial equation con¬ 
verges and agrees with the exact pressure from the energy 
equation, thus confirming the accuracy of our numerical 
results. These checks involve integrals of g 2 {r) and S{k) 
that are only slowly converging. Therefore, passing them 
demonstrates that our results have very high precision. 
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in Fig. 12. In the same dimension, as x increases, the 
distribution of Voronoi cell volumes narrows. This is 
expected because the system becomes more ordered as 
X increases. For the same y, the distribution also nar¬ 
rows as the dimension increases, consistent with theo¬ 
retical results that at fixed y, the nearest-neighbor dis¬ 
tance distribution narrows as dimension increases [12]. In 
Fig. 12, we additionally show the Voronoi cell-volume dis¬ 
tribution of saturated random sequential addition (RSA) 
[65-67] packings, the sphere packings generated by ran¬ 
domly and sequentially placing spheres into a large vol¬ 
ume subject to the nonoverlap constraint until no addi¬ 
tional spheres can be placed. Saturated RSA packings 
are neither stealthy nor hyperuniform [66, 67]. However, 
the Voronoi cell-volume distributions of saturated RSA 
packings look similar to that of the entropically favored 
stealthy ground states. This is not unexpected because 
Voronoi cell statistics are local characteristics, and hence 
are not sensitive to the stealthiness, which is a large-scale 
property. 


FIG. 11. (Color online) Structure factor and pair correlation 
fnnction for d = 2 for 0.33 < y < 0.46, as obtained from 
simulations. 


For smaller y values, the maximum of the structure 
factor is at the constraint cutoff k = K'^. However, for 
higher y values, the maximum of S{k) is no longer at 
k = 1'^. To probe this transition we have calculated the 
structure factor in two dimensions for 0.33 < y < 0.46. 
The results are shown in Fig. 11. As y increases, the 
peak at fc = 1+ gradually decreases its height, while the 
subsequent peak gradually grows and engulfs the first 
peak. 

Besides pair statistics, other widely used characteri¬ 
zation of point patterns include certain statistics of the 
Voronoi cells [14, 60-62]. A Voronoi cell is the region 
consisting of all of the points closer to a specific parti¬ 
cle than to any other. We have computed the Voronoi 
tessellation of the entropically favored stealthy ground 
states using the dD Convex Hulls and Delaunay Trian¬ 
gulations package [63] of the Computational Geometry 
Algorithms Library [64]. Since the number density of 
the stealthy ground states depends on the dimension and 
y, we rescaled each configuration to unity density for 
comparison of the Voronoi cell volumes. The probability 
distribution function p{vc) of the Voronoi cell volumes 
(where Vc is the volume of a Voronoi cell) are shown 
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FIG. 12. (Color online) Voronoi cell-volume distribution for 
1 < d < 3 for 0.05 < X ^ 0.25. For the same dimension, 
the Voronoi cell-volume distribution becomes narrower when 
X increases. For the same Xi th® Voronoi cell-volume distri¬ 
bution also becomes narrower when dimension increases. We 
also present Voronoi cell-volume distributions of RSA pack¬ 
ings at saturation here. 


One interesting phenomenon is that as x increases and 
approaches 1/2, systems that are not sufficiently large 
can become crystalline. In Fig. 13, we show two snap¬ 
shots of MD simulations at x = 0.48. The smaller config¬ 
uration is crystalline. However, systems that are 4 times 
larger remain disordered at the same x s^nd temperature. 
Therefore, this strongly indicates that crystallization is a 
finite-size effect for x tending to 1/2 from below. 





FIG. 13. (Color online) (a) Low-temperature MD snapshot 
of a 126-particle system at y = 0.48; the ground-state con¬ 
figuration is crystalline, (b) MD snapshot of a 504-particle 
system at the same Tg and x; the system does not crystallize 
and is indeed disordered without any Bragg peaks. 


B. X ^ 0-5 region 

As explained in Sec. II, we perform MD-based sim¬ 
ulated annealing with Monte Carlo moves of the simu¬ 
lation box for X > 0-5, since this method works bet¬ 
ter with rough potential energy surface and can mitigate 
the finite-size effect. We performed this simulation at 
X = 0.55, X = 0.73, and x = 0.81 in two dimensions. 
The results are shown in Fig. 14. The resulting con¬ 
figuration is always triangular lattice. Even though the 
ground-state manifold in this x regime contains aperiodic 
“wavy” phases discovered previously [14] [but which are 
called “stacked-slider” phases in the sequel to this pa¬ 
per [48], since they are aperiodic configurations with a 
high degree of order in which rows (in two dimensions) 
or planes (in three dimensions) of particles can slide past 
each other] as well as crystals other than the triangular 
lattice, the entropically favored ground state is always a 
triangular lattice. This means that the triangular lattice 
has a higher entropy than stacked-slider phases, although 
the latter appear to be more disordered [68]. 

Although we will show analytically that crystals are 
more entropically favored than stacked-slider phases in 
the upcoming paper of this series, we still need simula¬ 
tion results to determine which crystal structure has the 
highest entropy. The results of MD-based simulated an¬ 
nealing with Monte Carlo moves of the simulation box 
suggest that triangular lattice has the highest entropy 
in two dimensions. It seems natural to apply the same 
technique to three dimensions to determine the entropi- 
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cally favored crystal structure. However, we were unable 
to crystallize the system in three dimensions. Even the 
longest cooling schedule that we tried resulted in stacked- 
slider phases. 



FIG. 14. (Color online) MD-based simulated annealing result 
at (a) X = 0.55, (b) y = 0.73, and (c) y = 0.81. The ending 
configuration is triangular lattice except for small deforma¬ 
tions in the x ~ 0.55 case. 

Another way to find the entropically favored crystal is 
to use Wang-Landau Monte Carlo to directly calculate 
the entropy of different crystal structures as a function 
of the potential energy. We have performed this simula¬ 
tion on two-dimensional triangular lattice, square lattice, 
and three-dimensional body-centered cubic (BCC) lat¬ 
tice, face-centered cubic (FCC) lattice, and simple cubic 
(SC) lattice. The results are shown in Figs. 15 and 16. In 
all cases the entropy decreases as the energy decreases. In 
two dimensions, the entropy of the square lattice clearly 
decreases faster than that of the triangular lattice at 
every % value, confirming that the triangular lattice is 
entropically favored over the square lattice in the zero- 
temperature limit. In three dimensions at x = 0.58, the 
entropy of the FCC lattice decreases more slowly than 
that of the BCC and SC lattice, suggesting that the en¬ 
tropically favored ground state in three dimensions at 
X = 0.58 is the FCC lattice. At higher x values, the scal¬ 
ing of the entropy of the FCC lattice and the BCC lattice 


become very close to each other, preventing us from de¬ 
termining the entropically favored ground state at these 
X values. 



FIG. 15. (Golor online) Microcanonical entropy as a function 
of energy §($) calculated from Wang-Landau Monte Garlo of 
triangular lattice and square lattice at various x’s- Here "Fo 
denotes the ground-state energy. 



FIG. 16. (Golor online) Microcanonical entropy as a function 
of energy §($) calculated from Wang-Landau Monte Garlo of 
BCG lattice, FCC lattice, and SC lattice at various x’s. A 
curve for SC lattice is not presented for x > 0.68 because the 
latter is not a ground state at such high x values. Here "Fo 
denotes the ground-state energy. 


V. CONCLUSIONS AND DISCUSSION 

The uncountably infinitely degenerate classical ground 
states of the stealthy potentials have been sampled previ¬ 
ously using energy minimizations. We demonstrate here 
that this way of sampling the ground states to produce 
ensembles of configurations introduces dependencies on 
the energy minimization algorithm and the initial con¬ 
figuration. Such artificial dependencies are avoided in 
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studying the canonical ensemble in the T —>■ 0 limit. We 
sample this ensemble by performing MD simulations at 
sufficiently low temperatures, periodically taking snap¬ 
shots, and minimizing the energy of the snapshots. 

The configurations in this ensemble become more or¬ 
dered as X increases and obey certain theoretical condi¬ 
tions on their pair statistics [12], similarly to previous 
energy minimization results. However, other properties 
of this ensemble are unique. First, our numerical results 
demonstrate that the pair statistics of this ensemble dis¬ 
plays no “clustering effect” [divergence of g 2 (r) as r —> 0] 
for any y value, and is independent of the functional form 
of the stealthy potential. Second, we numerically verify 
the theoretical ansatz [12] that for sufficiently small x 
stealthy disordered ground states behave like “pseudo” 
disordered equilibrium hard-sphere systems in Fourier 
space, i.e., S{k) has the same functional form as the pair 
correlation function for equilibrium hard spheres for suf¬ 
ficiently small densities. Third, when y is above the crit¬ 
ical value of 0.5, our results strongly indicate that crystal 
structures are entropically favored in both two and three 
dimensions in the infinite-volume limit. Our numerical 
evidence suggests that the entropically favored crystal 
in two dimensions is the triangular lattice. However, 
we could not determine the entropically favored crystal 
structure in three dimensions. For finite systems, the 
disordered-to-crystal phase transition can happen at a 
slightly lower y. A theoretical explanation of this phe¬ 
nomenon remains an open problem. 

Besides ground states of stealthy potentials, other dis¬ 
ordered degenerate ground states of many-particle sys¬ 
tems have been studied using energy minimizations. 
Specifically, previous researchers have constrained the 
structure factor to have some targeted functional form 
other than zero for prescribed wave vectors [17, 18, 21]. 
Finding the configurations corresponding to such tar¬ 
geted structure factors amounts to finding the ground 
states of two-, three- and four-body potentials, in con¬ 
trast to the two-body stealthy potential studied in the 
present paper. This situation is the most general appli¬ 
cation of the collective-coordinate approach. It will be 
interesting to study the resulting pair statistics of the 
ground states for these more general interactions in the 
zero-temperature limit of the canonical ensemble. 

The collective-coordinate approach is an independent 
and fruitful addition to the basic statistical mechan¬ 
ics problem of connecting local interactions to macro¬ 
scopic observables. One important feature of collective- 
coordinate interactions is that it has uncountably in¬ 
finitely degenerate classical ground states [12]. In the 
case of isotropic pair interactions, the only other sys¬ 
tem that we know with this feature is the hard-sphere 
system. However, there are two important differences 
between hard-sphere systems and collective-coordinate 
ground states. First, while the dimensionality of the 
configuration space of equilibrium hard-sphere systems 
consisting of N particles within a periodic box is fixed 
[simply determined by the nontrivial number of degrees 


of freedom, 1)], the dimensionality of the collective- 

coordinate ground-state configuration space decreases as 
X increases and, on a per particle basis, eventually van¬ 
ishes [12]. The decreased dimensionality of the ground- 
state configuration space creates challenges for accurate 
sampling of the entropically favored ground states us¬ 
ing numerical simulations and hence the development of 
better sampling methods is a fertile ground for future 
research. 

Second, while the probability measure of the equi¬ 
librium hard-sphere system is uniform over its entire 
ground-state manifold, that of the stealthy ground states 
is not uniform. To illustrate this point, imagine a one¬ 
dimensional energy landscape that has a double-well po¬ 
tential behavior in a portion of the configuration space, 
as shown in Fig. 17. Each minimum represents a degen¬ 
erate ground state (as we find with stealthy potentials) 
and therefore the well depths of the minima are the same. 
Let us now consider harmonic approximations of the two 
wells in the vicinity of xi and X 2 , respectively, 

Vi{x) = ai{x - xif, 

and 

V2{x) = a2{x - X2f, 

where x is the configurational coordinate. At very low 
temperature, to a good approximation, the system can 
only visit the part of the configuration space with energy 
less than e, and £ —>■ 0 as T —)• 0. Solving Vi{x) < £, 
where i = 1,2, one finds the feasible region of configura¬ 
tion space associated with both wells: 

xi — \fejai < X < xi + \feja\, 

and 

a ^2 — \/ £/a 2 <x <X 2 + \fefa 2 - 

When ai 7 ^ 02 , we see that the feasible regions associated 
with the two potential wells have different ranges. There¬ 
fore, the weights associated with the two minima, i.e., the 
relative probabilities for finding the system in the vicinity 
of those minima, will also differ. Similarly, in the stealthy 
multidimensional configuration space that we are study¬ 
ing, the magnitude of the eigenvalues of the Hessian ma¬ 
trix will determine the relative weights. Therefore, the 
probability measure of the stealthy ground states is not 
uniform over the ground-state manifold, unlike the de¬ 
generate ground states of classical hard spheres. Our 
low-temperature MD simulations sample ground states 
with this nonuniform probability measure. It would be 
useful to devise theories to estimate the weights of dif¬ 
ferent portions of the ground-state manifold. However, a 
feature that complicates the problem is that the Hessian 
matrix has zero eigenvalues. In the associated directions 
of the eigenvectors of the configuration space, the energy 
scales more slowly than quadratically (harmonically) but 
we do not know the specific form. 
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Configurational Coordinate, x 


FIG. 17. (Color online) A model one-dimensional energy land¬ 
scape with two wells located at x\ and X 2 of the same depth 
but different curvatures. The “feasible regions,” i.e., regions 
where V{x) < e, is marked by red dashed lines. 

This paper, which investigates the entropically favored 
ground states, is the first of a two-paper series. In the 
second paper, we will study aspects of the ground-state 
manifold with an emphasis on configurations that are not 
entropically favored for x above 1 /2 (the ordered regime). 
In particular, we will more fully investigate the nature 
of so-called “wavy” crystals or “stacked-slider” phases, 
discovered in Ref. 14. Using an analytical description of 
such states, we will demonstrate that they are part of 
the ground state but are not entropically favored. Our 
analytical model will also demonstrate that stacked-slider 
phases exist in three and higher dimensions. 
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cells in two dimensions. Similarly, in three dimensions, 
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lattice is the unique ground state at Xmax- 

Appendix B: Local Gradient Descent Algorithm 

Most optimization algorithms are designed for effi¬ 
ciency. They use complex rules to determine the direc¬ 
tion of the next step and take as large steps as possible. 
These features make their path less obvious. To mini¬ 
mize energy in the path following the gradient vector, we 
designed a “local gradient descent algorithm” with the 
following steps: 

1. Start from an initial guess, x, and find the function 
value /(x) and derivative f'{x). 

2. Start from a relatively large (10“^ times the simula¬ 
tion box side length) step size, s, and calculate the 
vector to the next step Ax = ~s | j/[x)| ■ Find the 
function value at the next step /(x -|- Ax). Calcu¬ 
late the change of function value A/ = /(x-|-Ax) — 
/(x). 

3. If we are following the path of steepest descent ac¬ 
curately, the change of the function value should 
be close to /'(x) • Ax. If the difference between 
A/ and /'(x) • Ax is less than 1%, we accept this 
move. Otherwise, we abort this move and half the 
step size s. 

4. Repeat the above steps until a minimum is found 
with enough precision. 
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Appendix A: Real-space potential in finite systems 

In the infinite-system-size limit, an isotropic 'O(k) cor¬ 
respond to an isotropic real-space pair potential u(r). 
However, for finite systems, the corresponding v(r) is 
anisotropic. To illustrate the finite-size effect, we com¬ 
pare the two-dimensional real-space potential v(r) in the 
infinite-system-size limit to corresponding potentials as¬ 
sociated with finite-sized fundamental cells of square and 
rhombic shapes of different volumes in Fig. 18. The real- 
space potential in the rhombic simulation box with a 60° 
interior angle is appreciably more isotropic than the real- 
space potential in a square simulation box. Therefore, in 


Appendix G: Number of particles of every system in 
Sec. IV 


TABLE I. The number of particles N of each systems shown 
in Figs. 9 and 10. 


X 

V for d = 1 

iV for d = 2 

A for d = 3 

0.05 

1001 

541 

261 

0.1 

501 

270 

131 

0.143 

351 

190 

92 

0.2 

251 

136 

66 

0.25 

201 

109 

53 

0.33 

151 

181 

191 


In this appendix we report the number of particles N in 
each system in Sec. IV. Both configurations in Fig. 6 con¬ 
sist of 51 particles. Configurations (a) and (b) in Fig. 7 
consist of 271 and 151 particles, respectively. Those in 
Fig. 8 consist of 131 and 161 particles, respectively. 
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FIG. 18. (Color online) A portion of the real-space potential v(r) around the origin for the stealthy potential (5) with K — 1 
and V(k) — 1. (a)-(f) Real-space potential in a periodic simulation box that is [(a), (c), and (e)] square or [(b), (d), and (f)] 
rhombic in shape; the latter has a 60° interior angle. The volumes of the simulation boxes, v_f, are [(a) and (b)] 100, [(c) and 
(d)[ 400, and [(e) and (f)) 1385. Panels (a)-(d) use unrealistically small simulation boxes and is intended to illustrate finite-size 
effect only, (g) The real-space potential in the infinite-system-size limit. All potentials are normalized by their respective values 
at the origin since scaling does not affect the ground state. Note that, starting from the center, the dark (red) region indicates 
the highest values of the potential, whereas towards the edge of the box, the dark (blue) region indicates the lowest values of 
the potential. 


The number of particles of each system in Figs. 9, 10, 
11, and 12 are shown in Tables I, II, and III, respectively. 


Each configuration in Figs. 14, 15, and 16 consist of 36, 
400, and 343 particles, respectively. 
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